Interval Truth Model

Visualizations of Priors in the ITM

Authors

Affiliation

Matthias Kloft

University of Marburg

Björn S. Siepe

University of Marburg

Daniel W. Heck

University of Marburg

1 Latent Appraisal

1.1 True Intervals: Informative Beta Prior on Marginal Locations and Widths

1.1.1 Sample From Priors

# priors on marginal locations and widths
df <- data.frame(Tr_wid_splx = rbeta(1e5, 1.2, 3))
# multiplicative shifting parameters
s <- rbeta(1e5,1,1)
# lower bounds
df$Tr_L_splx <- 
  (1 - df$Tr_wid_splx) * s
# transform from simplex to bivariate normal
Tr_bvn <- simplex_to_bvn(
  cbind(
    df$Tr_L_splx,
    df$Tr_wid_splx,
    df$Tr_L_splx + df$Tr_wid_splx))

df <- cbind(df, Tr_loc = Tr_bvn[, 1], Tr_wid = Tr_bvn[, 2])

# transform back to simplex
Tr_splx <- bvn_to_simplex(
  cbind(df$Tr_loc, df$Tr_wid))
df$Tr_loc_splx <- Tr_splx[,1] + .5 * Tr_splx[,2]
df$Tr_wid_splx <- Tr_splx[,2]

1.1.2 Marginal Distributions

df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid),
               fill = "purple",
               alpha = .3,
               binwidth = .005, ) +
  scale_x_continuous(
    limits = c(NA, NA),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

1.1.3 Joint Distribution

1.1.4 Ternary

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc_splx, y = Tr_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

1.1.4.1 Unbounded Scale

df %>%
  ggplot() +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 1,
    alpha = .1,
    shape = 16
  ) +
  theme_pubr()

1.2 Lambda (Item Discernibility)

1.2.1 sigma_lambda

df$sigma_lambda <- exp(rnorm(1e5, log(.5), .5))
#median(lambda$sigma_lambda)

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_lambda),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.2.2 lambda

df$lambda_loc <- exp(rnorm(1e5, 0, df$sigma_lambda))
df$lambda_wid <- exp(rnorm(1e5, 0, df$sigma_lambda))

df %>%
  ggplot() +
  geom_histogram(aes(x = lambda_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3 E (Person Proficiency)

1.3.1 mu_E

Then narrow Student t prior ensures that the population proficiency is regularized to 1.

df$mu_E <- rnorm(1e5)

df %>%
  ggplot() +
  geom_histogram(aes(x = exp(mu_E)),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.2 sigma_E

df$sigma_E <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_E),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.3 E

df$E_loc <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))
df$E_wid <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))

df %>%
  ggplot() +
  geom_histogram(aes(x = E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.4 Lambda and E joint

df %>%
  ggplot() +
  geom_histogram(aes(x = 1 / lambda_loc / E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "SD of Latent Appraisal", y = "Density") +
  theme_pubr()

1.4 Omega (Residual Correlation)

beta_param <- 2
df$omega <- rbeta(1e5, shape1 = beta_param, shape2 = beta_param)*2-1

df %>%
  ggplot() +
  geom_histogram(aes(x = omega),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-1, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.5 Latent Appraisals

1.5.1 Marginal Appraisals

A <- rbvnorm(1e5, 
             mean1 = df$Tr_loc, 
             mean2 = df$Tr_wid,
             sd1 = 1 / df$lambda_loc / df$E_loc,
             sd2 = 1 / df$lambda_wid / df$E_wid,
             cor = df$omega)

df$A_loc <- rnorm(1e5, mean = df$Tr_loc, sd = 1 / df$lambda_loc / df$E_loc)
df$A_wid <- rnorm(1e5, mean = df$Tr_wid, sd = 1 / df$lambda_wid / df$E_wid)

A_splx <- bvn_to_simplex(df[,c("A_loc", "A_wid")])
df$A_loc_splx <- A_splx[,1] + .5 * A_splx[,2]
df$A_wid_splx <- A_splx[,2]

df %>%
  ggplot() +
  geom_density(aes(A_loc_splx),
               fill = "lightblue",
               alpha = .3,
               shape = 16) +
  geom_density(aes(A_wid_splx),
               fill = "purple",
               alpha = .3,
               shape = 16) +
  scale_x_continuous(limits = c(0, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Appraisal Locations / Widths", y = "Density") +
  theme_pubr()

1.6 Joint Appraisals

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = A_loc_splx, y = A_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2 Biases

2.1 Shifting Bias

2.1.1 sigma_b

df$sigma_b <- exp(rnorm(1e5, log(.5), 1))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_b),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.1.2 b

df$b_loc <- rnorm(1e5, mean = 0, sd = df$sigma_b)
df$b_wid <- rnorm(1e5, mean = 0, sd = df$sigma_b)

df %>%
  ggplot() +
  geom_histogram(aes(x = b_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-5, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

b_splx <- bvn_to_simplex(df[,c("b_loc", "b_wid")])
df$b_loc_splx <- b_splx[,1] + .5 * b_splx[,2]
df$b_wid_splx <- b_splx[,2]

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = b_loc_splx, y = b_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2.2 a (Scaling Bias)

2.2.1 sigma_a

df$sigma_a <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_a),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.2.2 a

df$a_loc <- exp(rnorm(1e5, mean = 0, sd = df$sigma_a))

df %>%
  ggplot() +
  geom_histogram(aes(x = a_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

3 Responses

3.1 Marginal Responses

df$Y_loc <- df$A_loc * df$a_loc + df$b_loc
df$Y_wid <- df$A_wid + df$b_wid

Y_splx <- bvn_to_simplex(df[,c("Y_loc", "Y_wid")])
df$Y_loc_splx <- Y_splx[,1] + .5 * Y_splx[,2]
df$Y_wid_splx <- Y_splx[,2]

df %>%
  ggplot() +
  geom_histogram(
    aes(Y_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Y_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

3.2 Joint Responses

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Y_loc_splx, y = Y_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

---

4 Recovery Check

4.1 Define Parameter Bounds for Variances

# compute a benchmark for the mean and SD of the parameters
mean_benchmark <- simplex_to_bvn(c(.4, .2, .4), type = "ilr")
sd_benchmark_loc <- simplex_to_bvn(c(.98, .01, .01), type = "ilr")
sd_benchmark_wid <- simplex_to_bvn(c(.495, .01, .495), type = "ilr")
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2]
# SD forTr_loc
sigma_Tr_loc <- sd_benchmark_loc[1] / 4
# SD Tr_wid
sigma_Tr_wid <- abs(sd_benchmark_wid[2] - mean_benchmark[2]) / 4

# SDs for other parameters
sigma_lambda_E_loc <- .3
sigma_lambda_E_wid <- .3
sigma_a_loc <- .3
sigma_b_loc <- sigma_Tr_loc / 3
sigma_b_wid <- sigma_Tr_wid / 3

4.2 Generate Data for One Replication

n_respondents <- 50
n_items <- 20

sim_data <-
  generate_itm_data_sim_study(
    n_respondents = n_respondents,
    n_items = n_items,
    mu_Tr_loc = mu_Tr_loc,
    mu_Tr_wid = mu_Tr_wid,
    sigma_Tr_loc = sigma_Tr_loc,
    sigma_Tr_wid = sigma_Tr_wid,
    sigma_lambda_E_loc = sigma_lambda_E_loc,
    sigma_lambda_E_wid = sigma_lambda_E_wid,
    sigma_a_loc = sigma_a_loc,
    sigma_b_loc = sigma_b_loc,
    sigma_b_wid = sigma_b_wid,
    omega_beta = 3
  )
responses <- sim_data$responses

4.2.1 Plot Data

samples <- gather_values(
  lower = responses$x_L,
  upper = responses$x_U,
  item_id = responses$jj,
  step_size = 0.01
)
samples %>% 
  ggplot(aes(x = samples)) +
  geom_density(fill = "purple", alpha = .5) +
  facet_wrap(~item_id) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.2 Plot True Intervals as Error Bars

data.frame(
  idx = 1:n_items,
  lower = sim_data$parameters$Tr_L,
  upper = sim_data$parameters$Tr_U
) %>%
  ggplot(aes(y = idx)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.3) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
  labs(x = "Latent Truth",
       y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.3 Plot True Intervals as Ternary Plot

data.frame(
  idx = 1:n_items,
  Tr_L = sim_data$parameters$Tr_L,
  Tr_U = sim_data$parameters$Tr_U) %>% 
  mutate(
  Tr_loc = Tr_L + (Tr_U - Tr_L) / 2,
  Tr_wid = Tr_U - Tr_L
) %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 4,
    alpha = .5,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

4.3 Fit Model

4.3.1 Stan Data

### Stan data declaration
I <- length(unique(responses$ii))
J <- length(unique(responses$jj))
N <- nrow(responses)
ii <- responses$ii
jj <- responses$jj
nn <- c(1:N)
Y_splx <- cbind(responses$x_splx_1, responses$x_splx_2, responses$x_splx_3) %>%
  as.matrix()

#Y_splx[is.na(Y_splx)]

# sanity check: all rows sum to 1
# check <- apply(Y_splx, 1, sum) %>% as.data.frame()
# table(check)

## stan data list
stan_data <- list(
  I = I,
  J = J,
  N = N,
  ii = ii,
  jj = jj,
  nn = nn,
  Y_splx = Y_splx,
  link = 1
)

4.3.2 Compile Model

# Choose model to fit
model_name <- "itm_simulation_v2_beta"
# Compile model
model <-
  cmdstanr::cmdstan_model(
    stan_file = here("src", "models", paste0(model_name, ".stan")),
    pedantic = TRUE,
    quiet = FALSE
  )

4.3.3 Run Sampler

# number of MCMC chains
n_chains <- 4
# Run sampler
fit_prior_check <- model$sample(
  data = stan_data,
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 500,
  adapt_delta = .99,
  init = .1
)

Running MCMC with 4 parallel chains...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 213.4 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 218.9 seconds.
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 222.3 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 263.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 229.5 seconds.
Total execution time: 263.4 seconds.

# save fit_prior_check
fit_prior_check$save_object(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))

# load  fit
fit_prior_check <- readRDS(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))

---

4.3.4 Get Estimates

estimates_summary <- fit_prior_check$summary()

4.3.5 Sampler Diagnostics

fit_prior_check$diagnostic_summary()

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.6690134 0.7145780 0.6690593 0.6188122

# sampler diagnostics
fit_prior_check$sampler_diagnostics(format = "df") %>% 
  psych::describe(quant = c(.05,.95),) %>%
  round(2) %>%  
  as.data.frame() %>% 
  dplyr::select(median, min, Q0.05, Q0.95,  max) %>% 
  .[-c(7:9),]

              median    min  Q0.05  Q0.95    max
treedepth__     7.00   6.00   6.00   7.00   7.00
divergent__     0.00   0.00   0.00   0.00   0.00
energy__      540.30 469.52 502.22 579.02 624.73
accept_stat__   0.99   0.91   0.97   1.00   1.00
stepsize__      0.04   0.04   0.04   0.05   0.05
n_leapfrog__  127.00  63.00 127.00 127.00 255.00

# convergence diagnostics
convergence_summary <- 
  fit_prior_check$draws(format = "df") %>%
  summarise_draws(.x = ., "rhat", "ess_bulk", "ess_tail") %>%
  remove_missing() %>%
  select(-variable) %>%
  psych::describe(., quant = c(.05, .95)) %>%
  as.data.frame() %>%
  select(mean, median, Q0.05, Q0.95, min, max)

convergence_summary %>% round(3)

             mean   median    Q0.05    Q0.95     min      max
rhat        1.002    1.002    1.000    1.005   0.999    1.010
ess_bulk 1956.513 1861.096  938.414 3327.645 417.098 4357.585
ess_tail 1467.217 1476.239 1161.061 1756.521 682.151 2000.306

4.3.6 Effective sample size (ESS) & Rhat Plots

# color scheme
color_scheme_set(scheme = "purple")

# Effective sample sizes
plot_neff <-
  mcmc_neff_hist(bayesplot::neff_ratio(fit_prior_check), binwidth = .01) +
  labs(title = "A") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  )
# Rhat
plot_rhat <-
  bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit_prior_check)) +
  labs(title = "B") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  ) +
  yaxis_text(on = TRUE)

# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)

4.4 Parameter Plots

# color scheme
color_scheme_set(scheme = "purple")

true_parameters <- sim_data$parameters

4.4.1 Person Competence: Location

plot_E_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_loc"), true = true_parameters$E_loc) +
  labs(
    subtitle = "Person Competence: Location",
  ) 
plot_E_loc

4.4.2 Person Competence: Width

plot_E_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_wid"), true = true_parameters$E_wid) +
  labs(subtitle = "Person Competence: Width")
plot_E_wid

4.4.3 Person Scaling Bias: Location

plot_a_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("a_loc"), true = true_parameters$a_loc) +
  labs(subtitle = "Person Scaling Bias: Location") +
  theme(
    text = element_text(
      family = "serif",
      size = 16,
      colour = 1
    ),
    plot.title = element_text(size = 20,
                              face = "bold"),
    plot.margin = unit(c(.1, .1, .1, .1), "cm"),
    axis.title = element_text(size = 16, color = 1),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.title.y = element_text(margin = margin(r = 5)),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x = element_text(margin = margin(t = 3))
  )
plot_a_loc

4.4.4 Person Shifting Bias: Location

plot_b_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_loc"), true = true_parameters$b_loc) +
  labs(subtitle = "Person Shifting Bias: Location") 
plot_b_loc

4.4.5 Person Shifting Bias: Width

plot_b_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_wid"), true = true_parameters$b_wid) +
  labs(subtitle = "Person Shifting Bias: Width")
plot_b_wid

4.4.6 Item Discernability: Location

plot_lambda_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_loc"),
                       true = true_parameters$lambda_loc)
plot_lambda_loc

4.4.7 Item Discernability: Width

# plot
plot_lambda_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_wid"),
                       true = true_parameters$lambda_wid) +
  labs(subtitle = "Item Discernability: Width")
plot_lambda_wid

4.4.8 Residual Correlation

plot_omega <-
  mcmc_recover_scatter(x = fit_prior_check$draws("omega"),
                       true = true_parameters$omega) +
  xlim(-1, 1) +
  ylim(-1, 1) +
  labs(subtitle = "Residual Correlation")
plot_omega

4.4.9 Latent Truth Intervals

# plot
plot_Tr_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_loc"),
                       true = true_parameters$Tr_loc) +
  labs(subtitle = "Item True Location")
plot_Tr_loc

# plot
plot_Tr_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_wid"),
                       true = true_parameters$Tr_wid) +
  labs(subtitle = "Item True Width")

plot_Tr_wid

# ITM estimates
latent_truth_est_ilr <- data.frame(
  idx = 1:J,
  type = "estimated",
  Tr_L_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[1:J],
  Tr_wid_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[(J + 1):(J * 2)]
) %>%
  mutate(Tr_loc_est = Tr_L_est + Tr_wid_est / 2, 
         Tr_U_est = Tr_L_est + Tr_wid_est)

# compute means in the unbounded space
means_ilr <-
  cbind(jj, simplex_to_bvn(Y_splx)) %>%
  as.data.frame() %>%
  dplyr::group_by(jj) %>%
  dplyr::summarise(across(everything(),mean, na.rm = TRUE)) %>%
  dplyr::ungroup() %>%
  dplyr::select(-jj) %>%
  bvn_to_simplex() %>%
  mutate(
    idx = 1:J,
    type = "simple_mean",
    Tr_L_sm = x_1,
    Tr_U_sm = 1 - x_3
  ) %>%
  select(-c(x_1, x_2, x_3))


# true parameters
latent_truth_true <- data.frame(
  idx = 1:J,
  type = "true",
  Tr_L_true = true_parameters$Tr_L,
  Tr_U_true = true_parameters$Tr_U
)

df_interval_plot_ilr <- full_join(latent_truth_est_ilr, latent_truth_true) %>%
  full_join(means_ilr) %>% 
  mutate(type = factor(type, levels = c("estimated", "true")))

4.4.9.1 Interval Plot

cols <- c("True" = "grey70","ITM (ILR)" = "red", "Mean" = "blue")

df_interval_plot_ilr %>%
  ggplot(aes(y = idx)) +
    geom_errorbarh(
    aes(xmin = Tr_L_true, xmax = Tr_U_true, col = "True"), 
    height = 0,
    linewidth = 5
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_sm, xmax = Tr_U_sm, col = "Mean"),
    height = 0,
    linewidth = 3,
    alpha = .4
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_est, xmax = Tr_U_est, col = "ITM (ILR)"),
    height = 0,
    linewidth = 2,
    alpha = .5
    ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
    scale_y_discrete(
    expand = expansion(add = 1)
  ) +
  scale_color_manual(values = cols) +
  labs(x = "Latent Truth", y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

5 Prior vs. Posterior

5.1 True Intervals

latent_truth <- data.frame(
  Tr_loc_splx_post = fit_prior_check$draws("Tr_loc_splx") %>%
    unlist() %>%
    as.vector(),
  Tr_wid_splx_post = fit_prior_check$draws("Tr_wid_splx") %>%
    unlist() %>%
    as.vector()
) %>%
  mutate(
    jj = factor(rep(1:J, each = nrow(.) / J)),
    Tr_wid_splx_prior = rbeta(nrow(.), 1.2, 3),
    Tr_loc_splx_prior = (1 - Tr_wid_splx_prior) * rbeta(nrow(.), 1, 1) + Tr_wid_splx_prior / 2
  )

Tr_bvn <- simplex_to_bvn(
  cbind(
    latent_truth$Tr_loc_splx_prior - .5 * latent_truth$Tr_wid_splx_prior,
    latent_truth$Tr_wid_splx_prior,
    1 - (
      latent_truth$Tr_loc_splx_prior + .5 * latent_truth$Tr_wid_splx_prior
    )
  )
)
latent_truth <- cbind(latent_truth,
                      Tr_loc_prior = Tr_bvn[, 1], Tr_wid_prior = Tr_bvn[, 2])


latent_truth %>%
  ggplot() +
  geom_abline(intercept = 0,
              slope = 2,
              alpha = .7) +
  geom_abline(intercept = 2,
              slope = -2,
              alpha = .7) +
  geom_density_2d_filled(
    aes(x = Tr_loc_splx_prior, 
        y = Tr_wid_splx_prior),
    binwidth = .05, 
    alpha = .3, ) +
  geom_point(
    aes(x = Tr_loc_splx_post, y = Tr_wid_splx_post, col = jj),
    size = .3,
    alpha = .3,
    shape = 16
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  labs(x = "Location", y = "Width") +
  # supress legends
  guides(col = FALSE, fill = FALSE) +
  theme_pubr()

5.2 Hyper-Priors

5.2.1 sigma_lambda: Location

colors <- c("prior" = "lightblue", "posterior" = "purple")

sigma_lambda_loc <- data.frame(posterior = fit_prior_check$draws("sigma_J[1]", format = "list") %>% unlist() %>%  as.vector())
sigma_lambda_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_loc %>%
  ggplot() +
  geom_density(aes(prior), fill = "lightblue", alpha = 1) +
  geom_density(aes(exp(posterior * .5 + log(.5))), fill = "purple", alpha = .3) +
  scale_x_continuous(limits = c(NA, 3), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.2 sigma_lambda: Width

sigma_lambda_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_J[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_lambda_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.3 sigma_E: Location

sigma_E_loc <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.4 sigma_E: Width

sigma_E_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.5 mu_E: Location

mu_E_loc <- data.frame(
  posterior = fit_prior_check$draws("mu_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_loc$prior <- rnorm(nrow(sigma_lambda_loc))

mu_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.6 mu_E: Width

mu_E_wid <- data.frame(
  posterior = fit_prior_check$draws("mu_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_wid$prior <-  rnorm(nrow(sigma_lambda_loc))

mu_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()